load surface_time;
S=length(x);
timehistint=linspace(timehist(1),timehist(end),length(timehist)); dt=timehistint(2)-timehistint(1);

lowcut=0;
surfacehisthf=surfacehist;
for s=1:S
    %surfacehist=interp1(timehist,surfacehist,timehistint); surfacehist(1)=surfacehist(2);
    tempsurface=interp1(timehist,surfacehist(s,:),timehistint); tempsurface(1)=tempsurface(2);
    temptemp=interp1(timehist,avtemphist(s,:),timehistint);
    avtemphist(s,:)=bandpass(temptemp,lowcut,inf,dt);
    surfacehist(s,:)=bandpass(tempsurface,lowcut,inf,dt);
    %surfacehist(s,:)=surfacehist(s,:)-surfacehist(s,1);
    surfacehisthf(s,:)=bandpass(tempsurface,1/100,inf,dt);
    
    for d=1:length(depthrange)
        temptemp=interp1(timehist,temphist(s,:,d),timehistint);
        %temptemp=temptemp-mean(temptemp);
        temptemp=bandpass(temptemp,lowcut,inf,dt);
        temphist(s,:,d)=temptemp;
    end
end


[TIMEHIST X]=meshgrid(timehistint,x);
figure(1); clf; surf(X,TIMEHIST,-surfacehist); xlabel('Distance/m'); ylabel('Time/Ma'); set(gca,'zdir','rev');
figure(2); clf; pcolor(X,TIMEHIST,-surfacehist); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar; 
%ax=axis; ax(3)=20; axis(ax);
print(gcf, '-djpeg100','-r100',['surface_time']);
figure(3); clf; pcolor(X,TIMEHIST,-surfacehisthf); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar; 
%ax=axis; ax(3)=20; axis(ax);
print(gcf, '-djpeg100','-r100',['surfacehf_time']);

for d=1:length(depthrange)
    figure(10+d); clf; pcolor(X,TIMEHIST,temphist(:,:,d)); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar; title(['Temperature at ' num2str(depthrange(d)-40e3) ' m']);
    %ax=axis; ax(3)=20; axis(ax);
    print(gcf, '-djpeg100','-r100',['temp_time' num2str(10+d)]);
    close;
end

d=20; lintemp=reshape(temphist(:,:,d),1,S*length(timehist)); linsurf=reshape(surfacehist(:,:),1,S*length(timehist)); linsurfhf=reshape(surfacehisthf(:,:),1,S*length(timehist));
figure(30); clf; plot(lintemp,linsurf,'.','markersize',0.1); hold on; %plot(lintemp,linsurfhf,'g.','markersize',0.1);
figure(31); clf; ax=plotyy(timehistint,surfacehisthf(round(S/2),:),timehistint,temphist(round(S/2),:,d)'); set(ax(1),'ydir','rev'); linkaxes(ax,'xy');

elasticdef=surfacehist*0; dx=x(2)-x(1);
for t=1:length(timehist)
    q=temphist(:,t,d)*3e-5*3300*300000;
    %q=avtemphist(:,t)*3e-5*3300;
    tempdef=flex1d(q,7.2e+10,0.17,50000,dx);
    %elasticdef(:,t)=flex1d(q,7.2e+10,0.17,50000,dx);
    elasticdef(:,t)=bandpass(tempdef,1/100,inf,dt);
end
figure(32); clf; pcolor(X,TIMEHIST,elasticdef); shading flat; xlabel('Distance/m'); ylabel('Time/Ma'); colorbar; title(['Deflection from temperature at' num2str(depthrange(d)-40e3) ' m']);
linelasticdef=reshape(elasticdef(:,:),1,S*length(timehist)); figure(30); plot(linelasticdef,linsurf,'y.','markersize',0.1); plot(linelasticdef,linsurfhf,'c.','markersize',0.1);


surfacehistflip=[surfacehist; flipud(surfacehist)];
correlation=zeros(1,S); weight=zeros(1,S);
for refindex=1:S
    for s=refindex:refindex+S-1
        correlation(s-refindex+1)=correlation(s-refindex+1)+(surfacehistflip(s,:)-mean(surfacehistflip(s,:)))*(surfacehistflip(refindex,:)-mean(surfacehistflip(refindex,:)))';
        weight(s-refindex+1)=weight(s-refindex+1)+1;
    end
end
correlation=correlation./weight;
correlation=correlation/correlation(1);
subplot(2,2,2); plot(x,correlation);

surfacehistflip=[surfacehisthf; flipud(surfacehisthf)];
correlation=zeros(1,S); weight=zeros(1,S);
for refindex=1:S
    for s=refindex:refindex+S-1
        correlation(s-refindex+1)=correlation(s-refindex+1)+(surfacehistflip(s,:)-mean(surfacehistflip(s,:)))*(surfacehistflip(refindex,:)-mean(surfacehistflip(refindex,:)))';
        weight(s-refindex+1)=weight(s-refindex+1)+1;
    end
end
correlation=correlation./weight;
correlation=correlation/correlation(1);
subplot(2,2,4); plot(x,correlation);